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ABSTRACT 

The power spectrum (PS) of mass density fluctuations, independent of "bias- 
ing" , is estimated from the Mark III catalog of peculiar velocities using Bayesian 
statistics. A parametric model is assumed for the PS, and the free parameters 
are determined by maximizing the probability of the model given the data. The 
method has been tested using detailed mock catalogs. It has been applied to 
generalized CDM models with and without COBE normalization. 

The robust result for all the models is a relatively high PS, with P(k)Q 1 ' 2 = 
(4.8±1.5) x 10 3 (/i _1 Mpc) 3 at k = 0.1 /iMpc" 1 . An extrapolation to smaller scales 
using the different CDM models yields o^fi ' 6 = 0.88 ± 0.15. The peak is weakly 
constrained to the range 0.02 < k < 0.06 /iMpc -1 . These results are consistent 
with a direct computation of the PS (Kolatt & Dekel 1996). When compared to 
galaxy- density surveys, the implied values for (3 (= Q°- 6 /b) are of order unity to 
within 25%. 

The parameters of the COBE-normalized, flat CDM model are confined by 
a 90% likelihood contour of the sort Vth^n v = 0.8 ± 0.2, where \i = 1.3 and 
v = 3.4, 2.0 for models with and without tensor fluctuations respectively. For open 
CDM the powers are \i = 0.95 and v = 1.4 (no tensor fluctuations). A T-shape 
model free of COBE normalization yields only a weak constraint: V = 0.4 ± 0.2. 



Subject headings: cosmology: theory — cosmology: observation — dark matter 
- galaxies: clustering — galaxies: distances and redshifts — large scale structure 
of universe 
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1. INTRODUCTION 

In the standard picture of cosmology, the structure on large scales originated from 
small-amplitude initial density fluctuations that were amplified by gravitational instability. 
These initial fluctuations are assumed to be a Gaussian random field, solely characterized 
by its power spectrum. On large scales, the fluctuations are linear even at late times, so 
that the power spectrum preserves its original shape. This makes it a very useful statistics 
for large-scale structure. 

The power spectra of galaxy density were derived for many different samples, in two 
angular dimensions or in three dimensions from redshift space. Unfortunately, these power 
spectra correspond to objects that are not necessarily unbiased tracers of the underlying 
mass distribution, and it is the mass distribution that is directly related to theory (e.g. 
Dekel & Rees 1987 for a review on "galaxy biasing"). Clear evidence for this bias is 
provided by the fact that galaxies of different types are observed to cluster differently 
(e.g. Dressier 1980). It would therefore be naive to assume that any of the galaxy power 
spectra directly reflects the mass PS. Furthermore, in estimates of the galaxy PS from 
redshift surveys, uncertainties also arise when correcting for redshift distortions (Kaiser 
1987, Zaroubi and Hoffman 1996). For these reasons, one wishes to measure the mass PS 
directly from dynamical data, bypassing the complex galaxy-biasing issues and the need to 
correct for redshift distortions. In principle, such dynamical information can be provided 
by peculiar velocities, by gravitational lensing effects, or by fluctuations in the cosmic 
microwave background (CMB). In particular, the accumulating catalogs of galaxy peculiar 
velocities enable a direct determination of the mass PS under the natural assumption that 
the galaxies are unbiased tracers of the large-scale, gravitationally-induced velocity field. 

The PS is computed here from the Mark III catalog of peculiar velocities (Willick et 
al. 1995 WI; 1996a WII; 1996b Will), which consists of more than 3000 galaxies. It was 
compiled from several different data sets of spiral and elliptical/SO galaxies with distances 
inferred by the forward Tully-Fisher and D n —a methods. These data were re-calibrated and 
self-consistently put together as a homogeneous catalog for velocity analysis. The catalog 
provides radial peculiar velocities and inferred distances with errors on the order of 17—21% 
of the distance per galaxy, sampled nonuniformly out to distances of ~ 80 /i _1 Mpc from 
the Local Group. 

The catalog exists in two versions that differ by the method of correction for Malmquist 
bias. In the standard "G" version, the galaxies are first heavily grouped into ~ 1200 objects 
ranging from isolated field galaxies to rich clusters. The grouping reduces the non-linear 
noise on small scales and the resulting Malmquist bias. Then the data are systematically 
corrected for the Malmquist bias (see Dekel 1994). As a reference, we also use the "S" 
version, where a more straightforward systematic correction is applied to the field galaxies 
as singles, at the expense of larger noise on small scales. 

These data allow a reasonable recovery of the dynamical fields with ~ 12/i -1 Mpc 
smoothing in a sphere of radius ~ 60/i _1 Mpc about the Local Group, extending to ~ 
80/i _1 Mpc in certain regions. The POTENT method (Bertschinger & Dekel 1989; Dekel, 
Bertschinger & Faber 1990; Dekel 1994, 1997) attempts a recovery of the underlying density 
field with fixed Gaussian smoothing within this volume. In an associated paper, Kolatt and 
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Dekel (1996, KD) have computed the mass PS from the smoothed density field recovered 
by POTENT from Mark III. The limitations of the data introduce severe systematic errors, 
that were modeled via Monte-Carlo mock catalogs and then used to correct the measured 
PS. Since the KD results naturally involve uncertainties, an independent estimate of the 
PS, using a very different method, is useful. 

Our purpose is to estimate the mass PS directly from the peculiar velocities of the 
Mark III catalog, by means of a likelihood analysis. The non-local nature of the peculiar 
velocities, i.e. being influenced by the mass distribution in a whole neighborhood, allows 
one to probe scales somewhat larger than those probed by the density field. For example, 
the effect of a bulk velocity across the entire volume is not evident if only the density field 
is considered. For a similar reason, the velocity field is expected to obey linear theory 
better than the density field smoothed on a comparable scale, and to closer resemble a 
Gaussian field. Our approach here does not involve any explicit window function, weighting 
or smoothing, nor does it require artificial binning of the PS. In addition, it automatically 
underweights noisy, unreliable data. 

The data analyzed here are especially suited for Bayesian analysis. The sparse and 
inhomogeneous sampling of a random Gaussian field with Gaussian errors yields a multi- 
variate Gaussian data set. The corresponding posterior probability distribution function 
(PDF) is a multivariate Gaussian that is completely determined by the assumed PS and 
the assumed covariance matrix of errors. Under these conditions one can write the joint 
PDF of the model PS and the underlying velocity or density field, and then simultaneously 
estimate the PS model parameters and recover the "Wiener filter" solution of the fields 
(Zaroubi et al. 1995). In an associated paper (Zaroubi, Hoffman & Dekel 1996), we present 
the high-resolution fields recovered from this same data set using the PS derived here. 

To apply our method, the simplifying assumptions that have to be made are (a) that 
the peculiar velocities are drawn from a Gaussian field, (b) that their correlation function 
can be derived from the density PS using linear theory, and (c) that the errors are Gaussian 
and accurately estimated. The need to assume a parametric functional form for the PS is 
also a limitation; one can try to achieve flexibility by using a large number of parameters 
and a variety of functional forms, but at the risk of making the likelihood analysis unstable 
in some cases (§5). 

The method is described in §2, where the relation between the PS and the velocity 
correlation functions is specified, and the likelihood algorithm for computing the PS is 
described. The method is tested using a mock catalog in §3. The resultant power spectra 
are presented in §4, as derived from the Mark III data alone, and for generalized CDM 
models imposing COBE normalization. The associated constraints on the cosmological 
parameters are analyzed. Our conclusions are summarized and discussed in §5. 

2. METHOD 

2.1. Velocity Correlations 

The computation of the matter power spectrum from the peculiar velocity data by 
means of likelihood analysis requires a relation between the velocity correlation function 
and the power spectrum. Define the two-point velocity correlation (3 x 3) tensor by the 
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average over all pairs of points and r,- that are separated by r = — r i? 

^V(r) = (^(r^fo)) , (1) 

where i> M (i"j) is the \i component of the peculiar velocity at r^. In linear theory, it can be ex- 
pressed in terms of two scalar functions of r = |r| (Gorski 1988), parallel and perpendicular 
to the separation r, 

^V( r ) = + - *±(r)]r li r v . (2) 

The spectral representation of these radial correlation functions is 

tt2-c2(q\ poo 

^ll( r ) = ^M/ o P(k)K U{ (kr)dk, (3) 

where K±(x) = j\{x)/x and K^(x) = jo—2ji(x)/x, with ji (x) the spherical Bessel function 
of order /. The cosmological O dependence enters as usual in linear theory via /(O) ~ O 6 , 
and H is the Hubble constant. A parametric functional form of P(k) thus translates to a 
parametric form of ^ 

2.2. Likelihood Analysis 

Let m be the vector of model parameters and d the vector of N data points. Then 
Bayes theorem states that the posterior probability density of a model given the data is 

W)-™. (4) 

The denominator is merely a normalization constant. The probability density of the model 
parameters, V(m), is unknown, and in the absence of any other information we assume 
it is uniform within a certain range. The conditional probability of the data given the 
model, V{d\m), is the likelihood function, £(d|m). The objective in this approach, which 
is finding the set of parameters that maximizes the probability of the model given the data, 
is thus equivalent to maximizing the likelihood of the data given the model (Kaiser 1988; 
see also Jaffe & Kaiser 1994 for a first application to the Lauer & Postman data). 

The Bayesian analysis measures the relative likelihood of different models. An absolute 
frequentist measure of goodness of fit could be provided by the Chi-square per degree of 
freedom, which we use as a check on the best parameters obtained by the likelihood 
analysis. 

Assuming that the velocities form a Gaussian random field, the two-point velocity 
correlation tensor^ fully characterizes the statistics of the velocity field. Define the radial- 
velocity correlation (N x N) matrix Uij by Uij = rforj, where i and j refer to the data 
points. Let the inferred radial peculiar velocity at be ttj, with the corresponding error 
also assumed to be a Gaussian random variable. The observed correlation matrix is then 
Uij = Uij + ef5ij : and the likelihood of the data points is 

C= [(27r)^det(^)]- 1/2 exp l~^U^ u A . (5) 
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Given that the correlation matrix, Uij, is symmetric and positive definite, we can use 
the Cholesky decomposition method (e.g. Press et al. 1992) for computing the likelihood 
function (Eq. 5). The significant contribution of the errors to the diagonal terms makes 
the matrix especially well suited for decomposition. The calculation for a given choice of 
parameters and N ~ 1200 data points takes a few minutes on a Dec- Alpha workstation 
(of SpecFP92 ~ 150). 

The likelihood function of Equation 5 is the posterior PDF of the parameters m. It 
is a x 2 distribution (with N degrees of freedom) with respect to the N data points, but 
it is not necessarily a x 2 distribution with respect to the parameters. Therefore, the task 
of assigning accurate confidence levels to the parameters requires elaborate integrations 
over the volume encompassed by the equal-likelihood surfaces in parameter space. In 
the present paper we limit ourselves to a rough estimate of confidence levels by crudely 
approximating — 21n C as a x 2 distribution in parameter space. 

An important feature of this method, which distinguishes it from other methods (e.g. 
KD), is that the formal likelihood errors include both the distance measurement errors and 
the cosmic scatter due to the finite discrete sampling. 

Note that the grouping in the Mark III catalog also serves as a mean of smoothing over 
nonlinear velocities. The PS on large-scales should not be much affected by this grouping 
because when group galaxies enter the likelihood analysis as individuals, they enter with 
low weights compared to the weight of a whole group. For un-grouped galaxies, the noise 
term efdij is larger (relative to the signal Uij) in the observed correlation matrix Uij of 
Equation 5. 

Note also that the quantity that can be derived from peculiar velocity data via the 
linear approximation is / 2 (fi) P(k), where P(k) is the mass density PS (see Eq. 3). 

3. TESTING THE METHOD 

Careful testing of the method with realistic mock catalogs is essential in view of the 
large distance errors, the sparse and non-uniform sampling, the bias-correction procedures, 
and the possible non-linear and non-Gaussian effects. 

The mock Mark III catalogs are described in Kolatt et al. (1996). They are based on 
simulations whose initial conditions were extracted from a reconstruction of the smoothed 
real-space density field from the IRAS 1.2Jy redshift survey, taken back into the linear 
regime. Small-scale perturbations were added by means of constrained random realizations. 
The system was then evolved forward in time using an N-body simulation assuming = 1, 
and stopped at two alternative times, when the rms density fluctuation in a top-hat sphere 
of radius 8 /i _1 Mpc reached cr 8 = 0.7 and later when cr 8 = 1.12. 

The "galaxies" in the simulation were identified via a linear biasing scheme (b=1.35), 
and they were divided into 'spirals' and 'ellipticals' according to Dressler's morphology- 
density relation. The galaxies were assigned TF quantities (internal velocities and absolute 
magnitudes) that were Gaussianly scattered about an assumed TF relation, and were then 
"observed" following the selection criteria of the actual data sets that compose the Mark 
III catalog. The mock catalogs were grouped and corrected for biases just like the real 
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data, producing analogous G and S mock catalogs. 
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Figure la: Contour map of ln-likelihood in the a — A plane for a mock catalog based on the parametric model 
of Equation 6 with 5=9.25 /i _1 Mpc. Contour spacing is A[ln£] = — 1. The best-fit point is marked. 

Figure lb: The true power spectrum of the simulation (heavy solid), compared with the best-fit solution (heavy 
dashed), and two power spectra whose parameters lie on the innermost closed contour of Figure fa. 



The true PS of the mass in the simulation is well approximated by the functional form 

with A = 4.68 and 12.28 x 10 4 ( /i^Mpc) 4 , B = 8.3 and 9.25^ _1 Mpc, and a = 3.2 and 
2.8, for the as =0.7 and 1.12 cases respectively. 

The likelihood analysis was applied to the mock catalogs using the parametric func- 
tional form of Equation 6 as a prior. Only two parameters were allowed to vary at a 
time while the third parameter was kept fixed. Figure la shows, for example, a contour 
map of ln£ (relative to the maximum-likelihood peak) for one of the as = 1.12, G mock 
catalogs, spanning the a — A plane with B = 9.25 /i -1 Mpc. The contours are separated 
by A[ln£] = —1. Maximum likelihood is obtained at A = A and a = 2.75 (compared to 
2.80). Assuming a % 2 distribution with two degrees of freedom, the 90% confidence limit 
of the likelihood around the best-fit parameters is at ln£ ps —2.3. We conclude that the 
recovery method works well. 

Figure lb shows the recovered PS in comparison with the true PS of the simulation. 
They almost coincide over the whole range of scales, showing slight deviations only on 
very small scales. To illustrate the level of uncertainty, we plot for comparison two other 
power spectra that were obtained with parameter pairs that lie on the innermost contour 
about the maximum in Figure la. It shows that the amplitude near the peak can be off 
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by about 25%, and that the recovery becomes more robust at moderately smaller scales. 
The success of the recovery is similar when the other pairs of parameters are allowed to 
vary, and also when allowing all three parameters to vary at the same time. 

Similar success was achieved when the method was applied to the S mock catalogs. 
In what follows we focus on the results from the G Mark III catalog, and refer to the S 
catalog as an indication for robustness. 



4. RESULTS 

4.1. TheT Model 

We first recover the PS from the velocity data alone, independent of COBE normal- 
ization. We use as a parametric prior the so-called V model (e.g. Efstathiou, Bond and 
White 1992), 



P(k) = AkT 2 (k), T(k) = (l + [ak/Y + (bk/Tf' 2 + (ck/T) 2 ]^ 



-l/u 



(7) 



with a = 6.4 /i _1 Mpc, b = 3.0 /i _1 Mpc, c = 1.7/i _1 Mpc and v = 1.13. The free parame- 
ters to be determined by the likelihood analysis are the normalization factor A and the V 
parameter. In the context of the CDM cosmological model, V has a specific cosmological 
interpretation, T = fi/i. Here, however, independently of CDM, Equation 7 serves as a 
generic function with logarithmic slopes n = 1 and —3 on large and small scales respec- 
tively, and with a turnover at some intermediate wavenumber that is determined by the 
single shape parameter V. 
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Figure 2a: Contour map of ln-likelihood for the T model. Contour spacing is A[ln£] = — 1. A in units of 
A„=2.0xl0 6 (/i~ 1 Mpc) 4 . 

Figure 2b: The most likely T- model power spectrum (solid), and five other models whose parameters lie on 
the innermost contour of Figure 2a. 
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Figure 2a shows the contour map of In C in the A — V plane. The maximum likelihood 
values are V = 0.5 ± 0.15 and A = 3.0 x 10 5 ( /i -1 Mpc) 4 . The corresponding value of 
crgfi ' 6 is 0.85 ±0.1. The error-bars quoted here (and throughout the paper) are the 90% 
confidence limits. Figure 2b shows the best-fit PS (solid). To illustrate the uncertainty 
in the PS we also show the power spectra of five other parameter pairs that lie on the 
innermost likelihood contour about the best fit (roughly the 65% confidence limit). We 
find the estimate of o^O ' 6 to be robust to the grouping of the data. On the other hand, 
the estimate of V without imposing COBE normalization is sensitive to the grouping; for 
the S catalog we obtain V = 0.3 ±0.1. 

4.2. COBE-Normalized CDM Models 

We now restrict our attention to the generalized family of CDM cosmological models, 
allowing variations in the cosmological parameters O, A and h, as well as the large-scale 
PS slope n and the contribution of tensor fluctuations. Furthermore, we now impose the 
normalization implied by the four-year COBE DMR data (Hinshaw et al. 1996) as an 
additional external constraint. The general form of the PS in these models is 

P(k) = A COBE {n, O, A) T 2 (0, Q B , h; k) k n , (8) 

where we adopt the CDM transfer function proposed by Sugiyama (1995, a slight modifi- 
cation of Bardeen et al. 1986): 

T(k) = ln ( 1 + 2 - 34 Q)) [! + 3. 899 + (16 . lg) 2 + {5Mq) s + (6-71?) 4] -V4 ^ (9fl) 

(96) 

(/i 50 = -£/o/50kms -1 Mpc -1 = 2h). The parameters are varied, two at a time, such 
that they span the range of currently popular CDM models, including Tilted-A CDM 
(flat: fl + A = 1, fl < 1, n < 1) and Tilted-Open CDM (A = 0, fl < 1, n < 1). 
We allow the possibility of nonzero tensor fluctuations, T/S = 7(1 — n), where the ratio 
is of quadrupole moments (C2) of tensor and scalar modes in the expansion of angular 
temperature fluctuations, (e.g. Turner 1993; Crittenden et al. 1993). In all cases, the 
baryonic density is assumed to be fi b = 0.024/i -2 (e.g. Tytler et al. 1996). 

The COBE normalization for each model has been calculated by various authors 
(Gorski et al. 1995; Sugiyama 1995; White & Bunn 1995), using different Boltzmann codes, 
different statistical analyses, and sometimes even different temperature maps. We have 
arbitrarily adopted Sugiyama's normalization as a backbone, and for models not studied 
by him we use the other results after matching them to Sugiyama's, using the models that 
they have investigated in common. 

In particular, the COBE normalization is modeled by Aqobe = Ai(Q)A2(n). For 
Tilted-A CDM models we use the fits: 

log Ai(fi) = 7.83 - 8.330 + 21.310 2 - 29.670 3 + 10.65O 4 + 15.420 5 - 

6.04O 6 - 13.970 7 + 8.610 8 , (10a) 



q = k 



Qh exp(-0 & - hl^Qb/Q) (/iMpc -1 ) 
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, , ( v / -2.78 + 2.78n T/5 = , inf . 

These fits are for h = 0.5, but the h dependence in the range of interest is weak, and we 
ignore it here. For Tilted-Open CDM model with T/S = the fit is: 

log Ai(fi) = 5.70 + 1.680 - 4.530 2 + 7.570 3 - 7.530 4 + 3.150 5 - 0.23O 6 , (11a) 

logA 2 (n) = -2.71 + 2.71n. (116) 
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h h 
Figure 3: Contour map of ln-likelihood in the h— CI plane for the ACDM (a) and OCDM (b) models with n=l. 
A[ln£] = — 1. Shown are the most likely values of il for two fixed values of h, and the error bars corresponding 
to 90% confidence level. 



4-2.1 Scale-Invariant Models 

Figure 3 shows the likelihood contour map, in the h — O plane, for the flat CDM 
(ACDM) and Open CDM (OCDM) families of models with n=l (four-year COBE nor- 
malization by Sugiyama). It is clear from the elongated contour pattern that O and h are 
not constrained very effectively independently of each other. It is a degenerate combina- 
tion of the two parameters that is being tightly determined by the elongated ridge of high 
likelihood. The constraints in the range plotted can be approximated (eye-ball fit) by the 
functional form 

fi/iso 1 - 3 = 0.83 + 0.09, ACDM, (12a) 
Oa 50 - 95 = 0.88 + 0.07, OCDM. (126) 

Also marked in Figure 3 are the most likely values of Q for fixed given values of h = 0.5 
and h = 0.75, and the corresponding one-dimensional error bars. The apparent weak 
preferences for high O (ACDM) or high h (OCDM) along the ridges of high likelihood is 
insignificant. In fact, it is not robust to changes between the G and S catalogs. 



9 





4-2.2 Tilted Models 

Figure 4 shows the likelihood in the n — O plane for the tilted- ACDM family of models. 
The computation is done for fixed values of h at h = 0.5 or 0.75, with or without tensor 
fluctuations. COBE normalization is by White and Bunn (1995). The elongated ridge of 
high likelihood can now be described by 

ft/iso 13 ™ 2 - = 0.83 ±0.12, T/S = 0, (13a) 
tthw^n 3 - 4 = 0.83 ±0.12, T/S = 7(1 - n) . (136) 

The trends along the ridges are, again, not robust to changing the data from G to S. 
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For a fixed Q, this relation can be understood qualitatively as follows: the normal- 
ization by COBE fixes the amplitude at small wavenumbers, k ~ 0.001, and the velocity 
data constrain the amplitude at k ~ 0.1. The wavenumber corresponding to the peak of 
the PS is roughly proportional to Vth. Therefore, if a good fit is obtained with certain 
values of h and n, a similarly good fit can be obtained with higher h and lower n, or 
vice versa. The presence of tensor fluctuations lowers the amplitude imposed by COBE at 
small wavenumbers, and thus weakens the requirement for a tilt in n. 




Peculiar Velocities + COBE 
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Figure 5: The same as Figure 4 for tilted Open CDM models without tensor component, A[ln £] = — !. 



Figure 5 shows the analogous likelihood results in the n — fi plane for the tilted- 
OCDM family of models (normalized by Sugiyama 1995). Computations were done for the 
two values of h as before, but this time only for the case of no tensor fluctuations. The 
results are qualitatively similar to the case of tilted-ACDM models. The tight constraint 
for tilted-OCDM is 

flh 50 - 95 n 1A = 0.88 ±0.09. (14) 
In this case, the trends along the ridge of high likelihood are weak. 



4-2.3 Model Dependence 

The most likely PS for the tilted-ACDM family of models with tensor fluctuations 
and h = 0.75, corresponding to £1 = 0.97 and n = 0.84, is plotted in Figure 6. The figure 
also shows the typical scatter in the PS about the most likely model, using n — O pairs 
that lie on the innermost likelihood contour (Fig. 4, lower-right panel). This scatter is 
somewhat smaller than in the V model (Fig. 2b), partly due to the additional constraint 
from COBE. Within this family of models, the constraints obtained are P(/c)0 12 = (4.2 ± 
0.8) x 10 3 (/i _1 Mpc) 3 at k = 0.1 /iMpc" 1 and a 8 O 6 = 0.88 ± 0.10, where the error-bars 
are the likelihood analysis 90% confidence limits. The scatter in these quantities for the 
other CDM families of models is roughly the same. 
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Figure 6: The PS of the most probable COBE-normalized CDM model (solid), and the scatter about it 
following parameter pairs that lie on the innermost likelihood contour (Fig. 4, lower-right panel). The COBE- 
free V model is also marked (heavy solid). The PS computed by KD from POTENT density of the same velocity 
data (independent of COBE or models), and their la measurement errors, are shown in three bins. 



Figure 7: The best-fit power spectra for various CDM models. 



How robust are the results to the choice of model within the CDM families? Table 1 
in Appendix A shows the features of the most likely models for each family of CDM variant 
discussed above, for the G catalog. Listed for each model is the corresponding value of 
cr 8 0°- 6 , the amplitude of the PS at k = 0.1 h Mpc -1 , the location of the peak and, most 
important, the functional fit to the high-likelihood ridge, in which the most likely model 
resides. 

These power spectra of the best-fit models are plotted in Figure 7. All the models 
agree to within ~20% for k > 0.1 /iMpc -1 , and they differ by up to 30 — 50% on larger 
scales. The amplitude of the PS at k = 0.1 h Mpc -1 across the models is constrained to 
be Pik)^ 1 - 2 = (4.8 ± 0.8) x 10 3 (/i^Mpc) 3 . The value of cr 8 0°- 6 is roughly the same for 
all these models and ranges between 0.87 and 0.94. 

Table 1 shows as well the ln£ value of each of these models (with the zero set arbi- 
trarily to the highest likelihood). The best-fit models are all of comparable likelihood. The 
most likely model, with a small margin, is the tilted-ACDM model with tensor fluctuations 
and h = 0.75: O = 0.97 and n = 0.84, the model plotted in Figure 6. This conclusion, 
however, is limited to the G catalog. For the S catalog, the preferred models are non-tilted 
(n = 1), with lower ft (~ 0.6 - 0.7) and h (~ 0.6) (see Table 2). 

Figure 6 also compares the PS corresponding to the most-likely COBE-normalized 
CDM model and the most likely V model that is independent of COBE. The agreement 
between the two is within ~ 10% on all scales. The agreement on the location of the 
peak in the PS reflects the fact that the peculiar velocity data themselves contain some 
meaningful information on scales as large as the size of the Mark III sample even without 
the constraint from COBE. 
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The PS computed by KD from the POTENT smoothed density field, that has been 
recovered from the same Mark III data independent of COBE or models, is also displayed 
in Figure 6. The results (for all the models tested here) agree within la of the measure- 
ment errors, and they agree particularly well near k = 0.1 h Mpc -1 , where the velocity 
data imposes the strongest constraints. The spectral slope m near k = 0.1 h Mpc -1 , as 
constrained by our method, is roughly m = —0.75 ±0.45. The KD spectrum appears to be 
steeper (m ~ —1-5), and the associated estimate of o^fi ' 6 by KD is therefore marginally 
lower, at ~ 0.7—0.8. In fact the KD slope is steeper than any of the CDM spectra discussed 
here (see Figure 7), and is roughly as steep as the PS predicted by the CHDM model, a 
7:3 mixture of cold and hot dark matter (KD). This may indicate that the CDM family, 
with the exception of CHDM, does not allow enough freedom for a perfect fit on small 
scales. Also, our present results may be less reliable on small scales because, other than 
grouping, we do not make here any correction for nonlinear effects. On the other hand, 
the result of the current paper is probably more reliable than KD on large scales, because 
the likelihood method uses all the velocity data including the large-scale flows, while the 
POTENT density field is insensitive to the bulk velocity. 




4-3. Robustness to Grouping in the Data 

To test robustness to the grouping and Malmquist correction in the Mark III data, 
we have performed a similar likelihood analysis to the S catalog. Table 2 in Appendix A 
is the equivalent of Table 1 and summarizes the results for the S catalog. The S likelihood 
analysis might be more susceptible to small-scale non-linear effects that are not properly 
dealt with in our method. Still, it is interesting to compare the results for the G and S 
catalogs in order to distinguish further between robust and non-robust conclusions. 
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Figure 8 shows, as an example, the results obtained from the G and S catalogs for 
the tilted-ACDM model with tensor fluctuations and h = 0.75. The resemblance between 
these results is typical of the other models as well. First, the high-likelihood ridges are 
similar in shape: the powers \i and v are almost the same in the G and S cases for all the 
CDM models we checked. Second, there is a slight difference in the constant defining the 
ridge: for the model plotted here the S result is O /iso 1 ' 25 n 3A = 0.71 ± 0.12 while the G 
ridge was Q h^Q 1 ' 3 n 3A = 0.85 ±0.12, preferring somewhat higher values for the parameters. 
A similar trend occurs for all the models (see Appendix A). 

The gradient of likelihood along the ridge is different for the G and S catalogs to 
the extent that it is sometimes of opposite sign. The G data slightly prefers a tilt in the 
spectrum and O, ~ 1, while the S data slightly favors lowering the value of O below unity 
and n ~ 1. This difference indicates that we should not assign high significance to the 
trends along the ridges of high likelihood and rather focus on the robust ridges themselves. 
We thus quote constraints on the degenerate combinations of O, h and n, and not on any 
of them separately. 

The likelihood analysis of the S data with the V model also yields some differences. The 
maximum likelihood value is V = 0.3 ±0.1 for the S catalog, compared to the V = 0.5 ±0.15 
for the G catalog. The shape parameter T (independent of COBE normalization) seems 
to be only poorly constrained by the data. Still, we can quote based on the velocity data 
alone that T > 0.2 at the 90% confidence level. 

The robust feature for the V model is the value of asfi 0,6 , which is stable at asfi ' 6 = 
0.85 ±0.1 for both the G and S catalogs. The lower T from the S data is accompanied by 
a higher amplitude A, such that the combination agQ ' 6 remains unchanged. 

Another robust feature, for all models, is the amplitude of the maximum-likelihood 
power spectrum on intermediate scales. The PS amplitude at k = 0.1 /iMpc" 1 for both S 
and G is constrained to be P(k)Q 1 - 2 = (4.8±1.5) x 10 3 (/i _1 Mpc) 3 . The location of the peak 
in the PS in the two cases is roughly the same, ranging between 0.025 and 0.04 h Mpc -1 
for the S catalog, and between 0.03 and 0.06 h Mpc -1 for the G catalog. 

5. CONCLUSION 

We have presented a Bayesian method for deriving the power spectrum of mass density 
fluctuations from the Mark III Catalog of Peculiar Velocities. The result is free of galaxy 
"biasing." The method extracts the maximum amount of useful information from the data. 
It is exact to first order, under the assumption of Gaussian fluctuations and Gaussian errors. 
Tests using realistic mock catalogs show that this approximation is adequate because of 
the large-scale coherence of the velocity field and because the large errors that dominate 
on small scales make the nonlinear effects contribute only weakly to the result. The formal 
likelihood errors include both measurement errors and cosmic scatter. Ignored are the 
errors in the COBE normalization, of order of a few percent. 

Our robust result for the whole family of models examined here as priors and for 
the different ways of handling the grouping of the data is that the mass PS amplitude at 
k = 0.1 h Mpc" 1 is P(/c)0 12 = (4.8 ±1.5) x 10 3 (/i _1 Mpc) 3 . The robust integral constraint 
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for the CDM family of models cx 8 fi - 6 = 0.88 ± 0.15. The errors quoted are crude: they 
reflect the typical 90% uncertainty for each of the best-fits within each family of models, 
combined with the typical scatter among these best-fit models and between the G and S 
treatments of the data. Similar results are obtained when using the velocity data alone 
and when the additional constraints from COBE are included. For the span of models 
checked, the PS peak is in the range 0.02 < k < 0.06 h Mpc -1 . 

This normalization of the PS is in pleasant agreement with the independent 'frequen- 
tist' computation of the PS via POTENT by KD, which yielded P(k)Q 1 - 2 = (4.6 ± 1.4) x 
10 3 (/i _1 Mpc) 3 at k = 0.1 h Mpc" 1 , and cr 8 fi - 6 « 0.7-0.8. Our new results differ at about 
the 2a level from the earlier estimate by Seljak & Bertschinger (1994) of cx 8 n - 6 = 1.3 ±0.3. 
They performed a likelihood analysis of the POTENT reconstructed density field based on 
the earlier Mark II sample. The main improvements since then are that the current anal- 
ysis includes five times denser sampling in a more extended volume, the systematic errors 
such as Malmquist bias are handled better, and a wider span of models and parameters is 
used in the likelihood analysis. It may be interesting to note that the current measurement 
is somewhat higher than the completely independent estimate of a similar quantity based 
on cluster abundances, a 8 O 56 ~ 0.57 ± 0.05 (White, Efstathiou, & Frenk 1993), but it is 
only about 2a away. 

The comparisons of the mass a 8 to the values observed for optical galaxies (~ 0.95) 
and for IRAS 1.2Jy galaxies (~ 0.6 — 0.7) indicate (3 values of order unity to within 25% 
for most galaxy types on these scales (see KD, Fig. 6 and Table 2, for more details). 

A T-shape model, free of COBE normalization, is constrained by the velocity data 
only weakly, to T in the range 0.2 — 0.6. The most likely value of T are somewhat higher 
than the canonical values of ~ 0.2 — 0.3 typically obtained from galaxy density surveys 
(e.g. Efstathiou et al. 1992, Peacock & Dodds 1994). 

Within the families of COBE-normalized CDM models, which we have restricted to the 
range O < 1 and n < 1, we have obtained constraints on combinations of the cosmo logical 
parameters of the sort O n v = 0.8 ±0.2, with \i = 1.3 for the flat models and with v = 
3.4 and 2.0 with and without tensor fluctuations respectively. For the open models without 
tensor fluctuations the corresponding powers are n = 0.95 and v = 1.4. The extended 
error quoted include the uncertainty due to measurement errors and finite sampling, the 
variation among the various CDM models, and the uncertainty in the best way to group 
the data for Malmquist correction. 

The most likely model for the G data is a tilted-ACDM model with tensor fluctuations, 
relatively high Hubble constant (h ~ 0.75), O near unity and n ~ 0.85. These results are 
consistent with the conclusion of White et al. (1995), who argue for tilted CDM based 
on several data sets including power spectra of galaxy density (Peacock & Dodds 1994), 
cluster correlations, pair- wise velocities and COBE's results. Their best fit is O = 1, 
h ~ 0.45, n = 0.9, with tensor fluctuations. This model is about la away from our best fit 
but the PS is quite similar. Our lower value of n compensates for the higher value of h. 

It is interesting to note that based on the velocity data and COBE normalization 
alone, the standard CDM model (0 = l,n = l,/i = 0.5) is less likely than its various 
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variants studied here (see Appendix A). A slight tilt in n or a small decrease of O can 
increase the likelihood significantly. However, this does not imply that the standard CDM 
model is necessarily ruled out as our analysis measures only relative likelihoods rather than 
absolute goodness of fit. 

In a subsequent paper (Zaroubi et al. 1996), our range of allowed power spectra from 
the peculiar velocities is translated to an angular power spectrum of CMB fluctuations and 
compared to recent subdegree observations. We find there that in order to also fit the height 
of the first CMB peak with a CDM spectrum, the tilt cannot be more pronounced than 
n ~ 0.9, the total energy density O + A should be close to unity, and the baryonic fraction 
should be large, ~ 0.1. A discussion of our results in the wider context of cosmological 
parameters can be found in Dekel, Burstein & White (1997) and Dekel (1997). 

It is worth recalling that the recovered power spectrum is highly sensitive to the 
assumed observational errors, because the noise has a systematic effect on the measured 
power spectrum. Effectively, the power spectrum of the errors has to be subtracted out. 
This can be done explicitly, as in KD, or implicitly, as in our current likelihood analysis. 
Our input error estimate is based on the careful error analysis of the Mark III data by 
Willick et al. (WI; WII; Will). If, for some reason, this is an underestimate of the errors, 
then the recovered power spectrum is an overestimate, and vice versa. 

Finally, a note of caution of technical nature about the method and its application to 
the current data under certain extreme conditions. If not enough constraints are imposed, 
the inversion of Uij (Eq. 5) may become dominated by the noise rather than the signal and 
might lead to a wrong answer. We encountered this problem when we tried to parameterize 
the PS with a multi-parameter function that did not enforce any upper bound on the power 
on large scales (bounds that are properly enforced when COBE constraints are used, for 
example). The likelihood analysis preferred in this case unphysically high power on large 
scales. This is probably, at least in part, a result of noise dominance, and is beyond the 
scope of this paper. An algorithm to detect and possibly eliminate this problem is discussed 
elsewhere (Zaroubi 1995). As long as we use models that are properly bound on large scales 
the results on small scales are robust. The V model, for example, does not suffer from this 
problem because its shape effectively enforces the appropriate bounds on large scales. The 
data strongly constrain the amplitude Po.i and the slope near k = 0.1 /iMpc -1 . A V model 
with the peak at a very small k and the same Pq.i would have required a steeper slope 
(m = —2 to —3) which is disfavored by the data. The V model thus enforces specific 
connections between properties that are strongly constrained by the data with properties 
that are poorly constrained. 
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APPENDIX A 

We present here the detailed likelihood analysis results for the various COBE normal- 
ized CDM families of models. Table 1 shows the results obtained for the G catalog. The 
functional fit of the high-likelihood ridge obtained in each case is given. The error-bar 
quoted is the formal 90% confidence limit of the degenerate parameter combination. The 
values of the cosmological parameters corresponding to the maximum-likelihood model 
within this ridge are listed as well. Parameters that were held fixed in the likelihood anal- 
ysis while varying the other parameters are denoted by "fixed" . For reference, the COBE 
normalized standard CDM model is also listed. The value of ln£ obtained for each of 
these models is given, with the zero set arbitrarily to correspond to the most likely model. 



TABLE 1 

Maximum-Likelihood Results of the G Catalog for CDM Models 









Po.i!J L! 














High Likelihood 




(h- 3 Mpc 3 ) 


(hMpc- 1 ) 










CDM Model 


Ridge 


(±0.10) 


(±1000) 


(±0.015) 


Q 


h 


n 


InC 


Standard 




1.2 


7700 


0.042 


1 (fixed) 


0.5(fixed) 


1 (fixed) 


-8.80 


A, n=l 


0/i 5 o L3 = 0.83 ±0.09 


0.92 


5600 


0.035 


0.97 


0.45 


1 (fixed) 


-2.34 


Open, n=l 


Qh 50 - 95 = 0.88 ±0.07 


0.87 


4600 


0.042 


0.51 


0.90 


1 (fixed) 


-0.45 


Tilted-A (T/S=0) 


fi/jso 1 - 3 ^ 2 ' 1 = 0.83± 0.12 


0.94 


5600 


0.037 


1.00 


0.5(fixed) 


0.92 


-1.97 




fi/i 5 o L3 n L9 = 0.83 ±0.12 


0.92 


4800 


0.044 


0.97 


0.75(fixed) 


0.72 


-0.78 


Tilted-A (T/S/0) 


fi/jso 1 - 3 ^ 3 ' 6 = 0.83 ± 0.11 


0.90 


5100 


0.039 


1.00 


0.5(fixed) 


0.95 


-1.19 




fi/i 50 L V' 4 = 0.85 ±0.11 


0.88 


4100 


0.053 


0.97 


0.75(fixed) 


0.84 


0.00 


Tilted-Open (T/S=0) 


nh 50 °- 95 n 1A = 0.88 ± 0.09 


0.91 


5300 


0.037 


1.00 


0.5(fixed) 


0.92 


-1.76 




fi/! 50 °' 9 V- 4 = 0.88 ±0.09 


0.91 


4700 


0.044 


0.97 


0.75(fixed) 


0.72 


0.70 



Several features of the power spectra of these best-fitting models are presented: the 
value of crgO ' 6 , the amplitude of the PS at k = 0.1 h Mpc -1 , Po.i^ 1 ' 2 , and the location 
of the peak in the PS, k pea k- The quoted error-bars in the table headings represent the 
typical 90% likelihood uncertainty in these quantities within each family of models. Due 
to the high-likelihood ridge in each family of models, the exact location of the maximum- 
likelihood model within the ridge is hardly significant, but the table demonstrates the 
robust features among these models (see discussion in §4.2.3). 
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Table 2 lists the results of the various CDM models for the S catalog. Similar high- 
likelihood ridges are listed, though there is a small systematic shift toward smaller values 
in the constant defining the ridge compared to the G data. The trends along the ridges 
and the likelihood ranking of the various models are different then in the G case but, again, 
they are only marginally significant. However the general features of the PS, such as the 
values of agO 6 and Pq.i^ 1 ' 2 , are similar (see §4.3). 



TABLE 2 

Maximum-Likelihood Results of the S Catalog for CDM Models 







<j 8 n - 6 




kpeak 












High Likelihood 




(h- 3 Mpc 3 ) 


(hMpc- 1 ) 










CDM Model 


Ridge 


(±0.10) 


(±1000) 


(±0.015) 


Q 


h 


?? 


InC 


Standard 




1.2 


7700 


0.042 


1 (fixed) 


0.5(fixed) 


1 (fixed) 


-40 


A, n=l 


Qho 1 - 25 = 0.69 ±0.07 


0.86 


5300 


0.030 


0.55 


0.60 


1 (fixed) 


-0.01 


Open, n=l 


O/i 50 a9 = 0.78 ±0.06 


0.83 


4800 


0.033 


0.90 


0.42 


1 (fixed) 


-0.24 


Tilted-A (T/S=0) 


fi/! 50 L2 V- 9 = 0.69 ± 0.09 


0.86 


5300 


0.031 


0.69 


0.5(fixed) 


1.00 


0.00 




fi/i50 L2 V- 9 = 0.70 ±0.09 


0.86 


5200 


0.033 


0.55 


0.75(fixed) 


0.88 


-0.20 


Tilted-A (T/S^O) 


fi/! 50 L25 n 3 - 4 = 0.69 ±0.12 


0.84 


5200 


0.031 


0.69 


0.5(fixed) 


1.00 


-0.02 




OAso 1 - 25 ^ 3 - 4 = 0.71 ± 0.12 


0.83 


4800 


0.032 


0.48 


0.75(fixed) 


0.96 


0.39 


Tilted-Open (T/S=0) 


fi/i 50 'V- 5 = 0.76 ±0.08 


0.84 


4800 


0.034 


0.76 


0.5(fixed) 


1.00 


-0.28 




fi/i 50 'V- 3 = 0.80 ±0.08 


0.83 


4500 


0.039 


0.55 


0.75(fixed) 


1.00 


-1.10 
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